/* Constructors for calculation */

#include "BaseFunction/basefunction.h"
/* constant --------------------------------------------------------------------------------------- */
#define SQR(x)     ((x)<0.0?-(x)*(x):(x)*(x))
/* observation code strings ----------------------------------------------------------------------- */
static string ObsCodes[] = {      

	""  ,"1C","1P","1W","1Y", "1M","1N","1S","1L","1E", /*  0- 9 */
	"1A","1B","1X","1Z","2C", "2D","2S","2L","2X","2P", /* 10-19 */
	"2W","2Y","2M","2N","5I", "5Q","5X","7I","7Q","7X", /* 20-29 */
	"6A","6B","6C","6X","6Z", "6S","6L","8L","8Q","8X", /* 30-39 */
	"2I","2Q","6I","6Q","3I", "3Q","3X","1I","1Q","5A"  /* 40-49 */
	"5B","5C","9A","9B","9C", "9X",""  ,""  ,""  ,""    /* 50-59 */
};
/* observations frequencies ----------------------------------------------------------------------- */
static unsigned char ObsFreqs[] = {
	/* 1:L1/E1, 2:L2/B1, 3:L5/E5a/L3, 4:L6/LEX/B3, 5:E5b/B2, 6:E5(a+b), 7:S */
	0, 1, 1, 1, 1,  1, 1, 1, 1, 1, /*  0- 9 */
	1, 1, 1, 1, 2,  2, 2, 2, 2, 2, /* 10-19 */
	2, 2, 2, 2, 3,  3, 3, 5, 5, 5, /* 20-29 */
	4, 4, 4, 4, 4,  4, 4, 6, 6, 6, /* 30-39 */
	2, 2, 4, 4, 3,  3, 3, 1, 1, 3, /* 40-49 */
	3, 3, 7, 7, 7,  7, 0, 0, 0, 0  /* 50-59 */
};
/* code priority table ---------------------------------------------------------------------------- */
static string CodePriors[7][MAXFREQ] = {  

	/* L1/E1      L2/B1        L5/E5a/L3 L6/LEX/B3 E5b/B2    E5(a+b)  S */
	{ "CPYWMNSL","PYWCMNDSLX","IQX"     ,""       ,""       ,""      ,"" }, /* GPS */
	{ "PC"      ,"PC"        ,"IQX"     ,""       ,""       ,""      ,"" }, /* GLO */
	{ "CABXZ"   ,""          ,"IQX"     ,"ABCXZ"  ,"IQX"    ,"IQX"   ,"" }, /* GAL */
	{ "CSLXZ"   ,"SLX"       ,"IQX"     ,"SLX"    ,""       ,""      ,"" }, /* QZS */
	{ "C"       ,""          ,"IQX"     ,""       ,""       ,""      ,"" }, /* SBS */
	{ "IQX"     ,"IQX"       ,"IQX"     ,"IQX"    ,"IQX"    ,""      ,"" }, /* BDS */
	{ ""        ,""          ,"ABCX"    ,""       ,""       ,""      ,"ABCX" }  /* IRN */
};

/* fatal callback function ------------------------------------------------------------------------ */
static fatalfunc_t *FatalFunc = NULL; 

/* crc tables generated by util/gencrc ------------------------------------------------------------ */
static const unsigned short Tbl_CRC16[] = {
	0x0000,0x1021,0x2042,0x3063,0x4084,0x50A5,0x60C6,0x70E7,
	0x8108,0x9129,0xA14A,0xB16B,0xC18C,0xD1AD,0xE1CE,0xF1EF,
	0x1231,0x0210,0x3273,0x2252,0x52B5,0x4294,0x72F7,0x62D6,
	0x9339,0x8318,0xB37B,0xA35A,0xD3BD,0xC39C,0xF3FF,0xE3DE,
	0x2462,0x3443,0x0420,0x1401,0x64E6,0x74C7,0x44A4,0x5485,
	0xA56A,0xB54B,0x8528,0x9509,0xE5EE,0xF5CF,0xC5AC,0xD58D,
	0x3653,0x2672,0x1611,0x0630,0x76D7,0x66F6,0x5695,0x46B4,
	0xB75B,0xA77A,0x9719,0x8738,0xF7DF,0xE7FE,0xD79D,0xC7BC,
	0x48C4,0x58E5,0x6886,0x78A7,0x0840,0x1861,0x2802,0x3823,
	0xC9CC,0xD9ED,0xE98E,0xF9AF,0x8948,0x9969,0xA90A,0xB92B,
	0x5AF5,0x4AD4,0x7AB7,0x6A96,0x1A71,0x0A50,0x3A33,0x2A12,
	0xDBFD,0xCBDC,0xFBBF,0xEB9E,0x9B79,0x8B58,0xBB3B,0xAB1A,
	0x6CA6,0x7C87,0x4CE4,0x5CC5,0x2C22,0x3C03,0x0C60,0x1C41,
	0xEDAE,0xFD8F,0xCDEC,0xDDCD,0xAD2A,0xBD0B,0x8D68,0x9D49,
	0x7E97,0x6EB6,0x5ED5,0x4EF4,0x3E13,0x2E32,0x1E51,0x0E70,
	0xFF9F,0xEFBE,0xDFDD,0xCFFC,0xBF1B,0xAF3A,0x9F59,0x8F78,
	0x9188,0x81A9,0xB1CA,0xA1EB,0xD10C,0xC12D,0xF14E,0xE16F,
	0x1080,0x00A1,0x30C2,0x20E3,0x5004,0x4025,0x7046,0x6067,
	0x83B9,0x9398,0xA3FB,0xB3DA,0xC33D,0xD31C,0xE37F,0xF35E,
	0x02B1,0x1290,0x22F3,0x32D2,0x4235,0x5214,0x6277,0x7256,
	0xB5EA,0xA5CB,0x95A8,0x8589,0xF56E,0xE54F,0xD52C,0xC50D,
	0x34E2,0x24C3,0x14A0,0x0481,0x7466,0x6447,0x5424,0x4405,
	0xA7DB,0xB7FA,0x8799,0x97B8,0xE75F,0xF77E,0xC71D,0xD73C,
	0x26D3,0x36F2,0x0691,0x16B0,0x6657,0x7676,0x4615,0x5634,
	0xD94C,0xC96D,0xF90E,0xE92F,0x99C8,0x89E9,0xB98A,0xA9AB,
	0x5844,0x4865,0x7806,0x6827,0x18C0,0x08E1,0x3882,0x28A3,
	0xCB7D,0xDB5C,0xEB3F,0xFB1E,0x8BF9,0x9BD8,0xABBB,0xBB9A,
	0x4A75,0x5A54,0x6A37,0x7A16,0x0AF1,0x1AD0,0x2AB3,0x3A92,
	0xFD2E,0xED0F,0xDD6C,0xCD4D,0xBDAA,0xAD8B,0x9DE8,0x8DC9,
	0x7C26,0x6C07,0x5C64,0x4C45,0x3CA2,0x2C83,0x1CE0,0x0CC1,
	0xEF1F,0xFF3E,0xCF5D,0xDF7C,0xAF9B,0xBFBA,0x8FD9,0x9FF8,
	0x6E17,0x7E36,0x4E55,0x5E74,0x2E93,0x3EB2,0x0ED1,0x1EF0
};
static const unsigned int Tbl_CRC24Q[] = {
	0x000000,0x864CFB,0x8AD50D,0x0C99F6,0x93E6E1,0x15AA1A,0x1933EC,0x9F7F17,
	0xA18139,0x27CDC2,0x2B5434,0xAD18CF,0x3267D8,0xB42B23,0xB8B2D5,0x3EFE2E,
	0xC54E89,0x430272,0x4F9B84,0xC9D77F,0x56A868,0xD0E493,0xDC7D65,0x5A319E,
	0x64CFB0,0xE2834B,0xEE1ABD,0x685646,0xF72951,0x7165AA,0x7DFC5C,0xFBB0A7,
	0x0CD1E9,0x8A9D12,0x8604E4,0x00481F,0x9F3708,0x197BF3,0x15E205,0x93AEFE,
	0xAD50D0,0x2B1C2B,0x2785DD,0xA1C926,0x3EB631,0xB8FACA,0xB4633C,0x322FC7,
	0xC99F60,0x4FD39B,0x434A6D,0xC50696,0x5A7981,0xDC357A,0xD0AC8C,0x56E077,
	0x681E59,0xEE52A2,0xE2CB54,0x6487AF,0xFBF8B8,0x7DB443,0x712DB5,0xF7614E,
	0x19A3D2,0x9FEF29,0x9376DF,0x153A24,0x8A4533,0x0C09C8,0x00903E,0x86DCC5,
	0xB822EB,0x3E6E10,0x32F7E6,0xB4BB1D,0x2BC40A,0xAD88F1,0xA11107,0x275DFC,
	0xDCED5B,0x5AA1A0,0x563856,0xD074AD,0x4F0BBA,0xC94741,0xC5DEB7,0x43924C,
	0x7D6C62,0xFB2099,0xF7B96F,0x71F594,0xEE8A83,0x68C678,0x645F8E,0xE21375,
	0x15723B,0x933EC0,0x9FA736,0x19EBCD,0x8694DA,0x00D821,0x0C41D7,0x8A0D2C,
	0xB4F302,0x32BFF9,0x3E260F,0xB86AF4,0x2715E3,0xA15918,0xADC0EE,0x2B8C15,
	0xD03CB2,0x567049,0x5AE9BF,0xDCA544,0x43DA53,0xC596A8,0xC90F5E,0x4F43A5,
	0x71BD8B,0xF7F170,0xFB6886,0x7D247D,0xE25B6A,0x641791,0x688E67,0xEEC29C,
	0x3347A4,0xB50B5F,0xB992A9,0x3FDE52,0xA0A145,0x26EDBE,0x2A7448,0xAC38B3,
	0x92C69D,0x148A66,0x181390,0x9E5F6B,0x01207C,0x876C87,0x8BF571,0x0DB98A,
	0xF6092D,0x7045D6,0x7CDC20,0xFA90DB,0x65EFCC,0xE3A337,0xEF3AC1,0x69763A,
	0x578814,0xD1C4EF,0xDD5D19,0x5B11E2,0xC46EF5,0x42220E,0x4EBBF8,0xC8F703,
	0x3F964D,0xB9DAB6,0xB54340,0x330FBB,0xAC70AC,0x2A3C57,0x26A5A1,0xA0E95A,
	0x9E1774,0x185B8F,0x14C279,0x928E82,0x0DF195,0x8BBD6E,0x872498,0x016863,
	0xFAD8C4,0x7C943F,0x700DC9,0xF64132,0x693E25,0xEF72DE,0xE3EB28,0x65A7D3,
	0x5B59FD,0xDD1506,0xD18CF0,0x57C00B,0xC8BF1C,0x4EF3E7,0x426A11,0xC426EA,
	0x2AE476,0xACA88D,0xA0317B,0x267D80,0xB90297,0x3F4E6C,0x33D79A,0xB59B61,
	0x8B654F,0x0D29B4,0x01B042,0x87FCB9,0x1883AE,0x9ECF55,0x9256A3,0x141A58,
	0xEFAAFF,0x69E604,0x657FF2,0xE33309,0x7C4C1E,0xFA00E5,0xF69913,0x70D5E8,
	0x4E2BC6,0xC8673D,0xC4FECB,0x42B230,0xDDCD27,0x5B81DC,0x57182A,0xD154D1,
	0x26359F,0xA07964,0xACE092,0x2AAC69,0xB5D37E,0x339F85,0x3F0673,0xB94A88,
	0x87B4A6,0x01F85D,0x0D61AB,0x8B2D50,0x145247,0x921EBC,0x9E874A,0x18CBB1,
	0xE37B16,0x6537ED,0x69AE1B,0xEFE2E0,0x709DF7,0xF6D10C,0xFA48FA,0x7C0401,
	0x42FA2F,0xC4B6D4,0xC82F22,0x4E63D9,0xD11CCE,0x575035,0x5BC9C3,0xDD8538
};

/* functions -------------------------------------------------------------------------------------- */
/* decoded binary data to ASCII data -------------------------------------------------------------- */
/* extract unsigned/signed bits ----------------------------------------------------------------------
* extract unsigned/signed bits from byte data
* args   : unsigned char *ChBuff  I  byte data
*          int           BitPos	  I  bit position from start of data (bits)
*          int           BitLen   I  bit length (bits) (len<=32)
* return : extracted unsigned/signed bits
*-------------------------------------------------------------------------------------------------- */
unsigned int getbitu(const unsigned char *ChBuff,int BitPos,int BitLen){
	unsigned int nBits=0;
	int i;
	for (i=BitPos; i<BitPos+BitLen; i++) nBits=(nBits<<1)+((ChBuff[i/8]>>(7-i%8))&1u);
	return nBits;
}
int getbits(const unsigned char *ChBuff,int BitPos,int BitLen)
{
	unsigned int nBits=getbitu(ChBuff,BitPos,BitLen);
	if (BitLen<=0||32<=BitLen||!(nBits&(1u<<(BitLen-1)))) return (int)nBits;
	return (int)(nBits|(~0u<<BitLen)); /* extend sign */
}
/* set unsigned/signed bits --------------------------------------------------------------------------
* set unsigned/signed bits to byte data
* args   : unsigned char *ChBuff  IO byte data
*          int           BitPos   I  bit position from start of data (bits)
*          int           BitLen   I  bit length (bits) (len<=32)
*         (unsigned) int ByteData I  unsigned/signed data
* return : none
*-------------------------------------------------------------------------------------------------- */
void setbitu(unsigned char *ChBuff,int BitPos,int BitLen,unsigned int ByteData){
	unsigned int nMask=1u<<(BitLen-1);
	int i;
	if (BitLen<=0||32<BitLen) return;
	for (i=BitPos; i<BitPos+BitLen; i++,nMask>>=1) {
		if (ByteData&nMask) ChBuff[i/8]|=1u<<(7-i%8); else ChBuff[i/8]&=~(1u<<(7-i%8));
	}
}
/* crc-24q parity ------------------------------------------------------------------------------------
* compute crc-24q parity for sbas, rtcm3
* args   : unsigned char *ChBuff  I   data
*          int           BitLen   I      data length (bytes)
* return : crc-24Q parity
* notes  : see reference [2] A.4.3.3 Parity
* ------------------------------------------------------------------------------------------------- */
unsigned int rtk_crc24q(const unsigned char *ChBuff,int BitLen){
	unsigned int nCrc=0;
	int i;
	for (i=0; i<BitLen; i++) nCrc=((nCrc<<8)&0xFFFFFF)^Tbl_CRC24Q[(nCrc>>16)^ChBuff[i]];
	return nCrc;
}
/* crc-16 parity -------------------------------------------------------------------------------------
* compute crc-16 parity for binex, nvs
* args   : unsigned char *ChBuff I data
*          int           BitLen  I      data length (bytes)
* return : crc-16 parity
* notes  : see reference [10] A.3.
*-------------------------------------------------------------------------------------------------- */
unsigned short rtk_crc16(const unsigned char *ChBuff,int BitLen)
{
	unsigned short nCrc=0;
	int i;

	for (i=0; i<BitLen; i++) {
		nCrc=(nCrc<<8)^Tbl_CRC16[((nCrc>>8)^ChBuff[i])&0xFF];
	}
	return nCrc;
}
/* decode navigation data word -----------------------------------------------------------------------
* check party and decode navigation data word
* args   : unsigned int NavWord  I navigation data word (2+30bit)
*                              (previous word D29*-30* + current word D1-30)
*          unsigned char *ChData O decoded navigation data without parity
*                              (8bitx3)
* return : status (1:ok,0:parity error)
* notes  : see reference [1] 20.3.5.2 user parity algorithm
*-------------------------------------------------------------------------------------------------- */
int decode_word(unsigned int NavWord,unsigned char *ChData)
{
	const unsigned int Hammings[]={
		0xBB1F3480,0x5D8F9A40,0xAEC7CD00,0x5763E680,0x6BB1F340,0x8B7A89C0
	};
	unsigned int nParity=0,w;
	int i;

	if (NavWord&0x40000000) NavWord^=0x3FFFFFC0;

	for (i=0; i<6; i++) {
		nParity<<=1;
		for (w=(NavWord&Hammings[i])>>6; w; w>>=1) nParity^=w&1;
	}
	if (nParity!=(NavWord&0x3F)) return 0;

	for (i=0; i<3; i++) ChData[i]=(unsigned char)(NavWord>>(22-i*8));
	return 1;
}
void matmul_pnt(const char *TraFlag,int SizeA,int SizeB,int SizeAB,double CoeAB,
	const double *MatA,const double *MatB,double CoeC,double *MatC){
	double LineAB;
	int i,j,x,f=TraFlag[0]=='N' ? (TraFlag[1]=='N' ? 1 : 2) : (TraFlag[1]=='N' ? 3 : 4);

	for (i=0; i<SizeA; i++) for (j=0; j<SizeB; j++) {
		LineAB=0.0;
		switch (f) {
		case 1: for (x=0; x<SizeAB; x++) LineAB+= MatA[i+x*SizeA]  *MatB[x+j*SizeAB]; break;
		case 2: for (x=0; x<SizeAB; x++) LineAB+= MatA[i+x*SizeA]  *MatB[j+x*SizeB]; break;
		case 3: for (x=0; x<SizeAB; x++) LineAB+= MatA[x+i*SizeAB] *MatB[x+j*SizeAB]; break;
		case 4: for (x=0; x<SizeAB; x++) LineAB+= MatA[x+i*SizeAB] *MatB[j+x*SizeB]; break;
		}
		if (CoeC==0.0) MatC[i+j*SizeA] = CoeAB*LineAB;
		else MatC[i+j*SizeA] = CoeAB*LineAB + CoeC*MatC[i+j*SizeA];
	}
}
void matmul_vec(const char *TraFlag,int SizeA,int SizeB,int SizeAB,double CoeAB,
	const vector<double> &MatA,const vector<double> &MatB,double CoeC,vector<double> &MatC){
	double LineAB;
	int i,j,x,f=TraFlag[0]=='N' ? (TraFlag[1]=='N' ? 1 : 2) : (TraFlag[1]=='N' ? 3 : 4);

	for (i=0; i<SizeA; i++) for (j=0; j<SizeB; j++) {
		LineAB=0.0;
		switch (f) {
		case 1: for (x=0; x<SizeAB; x++) LineAB+= MatA[i+x*SizeA]  *MatB[x+j*SizeAB]; break;
		case 2: for (x=0; x<SizeAB; x++) LineAB+= MatA[i+x*SizeA]  *MatB[j+x*SizeB]; break;
		case 3: for (x=0; x<SizeAB; x++) LineAB+= MatA[x+i*SizeAB] *MatB[x+j*SizeAB]; break;
		case 4: for (x=0; x<SizeAB; x++) LineAB+= MatA[x+i*SizeAB] *MatB[j+x*SizeB]; break;
		}
		if (CoeC==0.0) MatC[i+j*SizeA] = CoeAB*LineAB;
		else MatC[i+j*SizeA] = CoeAB*LineAB + CoeC*MatC[i+j*SizeA];
	}
}
void matmul_vec(const char *TraFlag,int SizeA,int SizeB,int SizeAB,double CoeAB,
	const vector<double> &MatA,const vector<double> &MatB,double CoeC,vector<double> &MatC,
	const int StartA, const int StartB, const int StartC){
	double LineAB;
	int i,j,x,f=TraFlag[0]=='N' ? (TraFlag[1]=='N' ? 1 : 2) : (TraFlag[1]=='N' ? 3 : 4);

	for (i=0; i<SizeA; i++) for (j=0; j<SizeB; j++) {
		LineAB=0.0;
		switch (f) {
		case 1: for (x=0; x<SizeAB; x++) 
			LineAB+= MatA[StartA+i+x*SizeA]  *MatB[StartB+x+j*SizeAB]; break;
		case 2: for (x=0; x<SizeAB; x++) 
			LineAB+= MatA[StartA+i+x*SizeA]  *MatB[StartB+j+x*SizeB]; break;
		case 3: for (x=0; x<SizeAB; x++) 
			LineAB+= MatA[StartA+x+i*SizeAB] *MatB[StartB+x+j*SizeAB]; break;
		case 4: for (x=0; x<SizeAB; x++) 
			LineAB+= MatA[StartA+x+i*SizeAB] *MatB[StartB+j+x*SizeB]; break;
		}
		if (CoeC==0.0) MatC[StartC+i+j*SizeA] = CoeAB*LineAB;
		else MatC[StartC+i+j*SizeA] = CoeAB*LineAB + CoeC*MatC[StartC+i+j*SizeA];
	}
}
/* (static) LU decomposition ------------------------------------------------------ */
int LUdcmp(vector<double> &MatSrc,int Size,vector<int> &index){
	double big,s,tmp;
	int i,imax=0,j,k;
	vector<double> vv(Size,0.0);

	for (i=0; i<Size; i++) {
		big=0.0; 
		for (j=0; j<Size; j++) if ((tmp=fabs(MatSrc[i+j*Size]))>big) big=tmp;
		if (big>0.0) vv[i]=1.0/big; 
		else { vv.clear(); return -1; }
	}
	for (j=0; j<Size; j++) {
		for (i=0; i<j; i++) {
			s=MatSrc[i+j*Size]; 
			for (k=0; k<i; k++) s-=MatSrc[i+k*Size]*MatSrc[k+j*Size]; 
			MatSrc[i+j*Size]=s;
		}
		big=0.0;
		for (i=j; i<Size; i++) {
			s=MatSrc[i+j*Size]; 
			for (k=0; k<j; k++) s-=MatSrc[i+k*Size]*MatSrc[k+j*Size]; 
			MatSrc[i+j*Size]=s;
			if ((tmp=vv[i]*fabs(s))>=big) { big=tmp; imax=i; }
		}
		if (j!=imax) {
			for (k=0; k<Size; k++) {
				tmp=MatSrc[imax+k*Size]; 
				MatSrc[imax+k*Size]=MatSrc[j+k*Size]; MatSrc[j+k*Size]=tmp;
			}
			vv[imax]=vv[j];
		}
		index[j]=imax;
		if (MatSrc[j+j*Size]==0.0) { vv.clear(); return -1; }
		if (j!=Size-1) {
			tmp=1.0/MatSrc[j+j*Size]; 
			for (i=j+1; i<Size; i++) MatSrc[i+j*Size]*=tmp;
		}
	}
	vv.clear();
	return 0;
}
/* LU back-substitution ----------------------------------------------------------- */
void LUbksb(vector<double> &MatB,int Size,vector<int> &index,vector<double> &MatA,int startA){
	double s;
	int i,ii=-1,ip,j;

	for (i=0; i<Size; i++) {
		ip=index[i]; s=MatA[startA+ip]; MatA[startA+ip]=MatA[startA+i];
		if (ii>=0) for (j=ii; j<i; j++) s-=MatB[i+j*Size]*MatA[startA+j];
		else if (s) ii=i;
		MatA[startA+i]=s;
	}
	for (i=Size-1; i>=0; i--) {
		s=MatA[startA+i];
		for (j=i+1; j<Size; j++) s-=MatB[i+j*Size]*MatA[startA+j];
		MatA[startA+i]=s/MatB[i+i*Size];
	}
}
/* inverse of matrix -------------------------------------------------------------- */
int matinv(vector<double> &MatSrc,int Size){
	vector<double> matBuf;
	vector<int> index(Size,1);

	/* initialize matBuf use MatSrc */
	matBuf.assign(MatSrc.begin(),MatSrc.end());

	if (LUdcmp(matBuf,Size,index)==-1) { matBuf.clear(); index.clear(); return -1; }
	for (int j=0; j<Size; j++){
		for (int i=0; i<Size; i++) MatSrc[i+j*Size]=0.0;
		MatSrc[j+j*Size]=1.0;
		LUbksb(matBuf,Size,index,MatSrc,j*Size);
	}

	matBuf.clear(); index.clear();
	return 0;
}

/* polynomial coefficients estimate ----------------------------------------------- */
int polyest(const vector<double> &A,const vector<double> &L,vector<double> &Coe,
	const int numL,const int numX,double &sigma){
	if (numL<numX) return -1;
	vector<double> AA(numX*numX,0.0),AL(numX,0.0),V(L.begin(),L.end());

	matmul_vec("NT",numX,numX,numL,1.0,A,A,0.0,AA);//AA
	matmul_vec("NT",numX,1,numL,1.0,A,L,0.0,AL);//AL
	if (matinv(AA,numX)==-1) return -1;
	matmul_vec("NN",numX,1,numX,1.0,AA,AL,0.0,Coe);//Coe

																	   /* sigma */
	matmul_vec("TN",numL,1,numX,1.0,A,Coe,-1.0,V);//V
	sigma=norm(V.begin(),numL)/sqrt(numL>numX ? numL-numX : 1);

	return 1;
}
/* solve linear equation ---------------------------------------------------------- */
int solve_line(const string TransFlag,const vector<double> &A,const vector<double> &Y,
	vector<double> &X,int Xnum,int SolNum) {
	vector<double> A_(A.begin(),A.end());
	
	string trans = TransFlag[0]=='N'? "NN" : "TN";

	if (matinv(A_,Xnum)==0) 
		matmul_vec(trans.c_str(),Xnum,SolNum,Xnum,1.0,A_,Y,0.0,X);
	else return -1;

	A_.clear();
	return 0;
}
/* get index -----------------------------------------------------------------*/
int getindex(double value, const double *range){
	if (range[2]==0.0) return 0;
	if (range[2]>0.0&&(value<range[0]||range[1]<value)) return -1;
	if (range[2]<0.0&&(value<range[1]||range[0]<value)) return -1;
	return (int)floor((value-range[0])/range[2]+0.5);
}
/* get number of items -------------------------------------------------------*/
int nitem(const double *range){
	return getindex(range[1],range)+1;
}

/* data format transfer functions ----------------------------------------------------------------- */
/* replace string ----------------------------------------------------------------- */
static int repstr(string &BasStr,const string TarStr,const string RepStr){
	int i=0;
	size_t strPos;

	while ((strPos=BasStr.find(TarStr))!=string::npos){
		BasStr.replace(strPos,TarStr.length(),RepStr);
		i=1;
	}
	return i;
}
/* replace keywords in file path -----------------------------------------------------
* replace keywords in file path with date, time, rover and base station id
* args   : string  SrcPath     I   file path (see below)
*          string  DstPath     O   file path in which keywords replaced (see below)
*          gtime_t GpsTime     I   time (gpst)  (time.time==0: not replaced)
*          string  RovID       I   rover id string        ("": not replaced)
*          string  BaseID      I   base station id string ("": not replaced)
* return : status (1:keywords replaced, 0:no valid keyword in the path,
*                  -1:no valid time)
* notes  : the following keywords in path are replaced by date, time and name
*              %Y -> yyyy : year (4 digits) (1900-2099)
*              %y -> yy   : year (2 digits) (00-99)
*              %m -> mm   : month           (01-12)
*              %d -> dd   : day of month    (01-31)
*              %h -> hh   : hours           (00-23)
*              %M -> mm   : minutes         (00-59)
*              %S -> ss   : seconds         (00-59)
*              %n -> ddd  : day of year     (001-366)
*              %W -> wwww : gps week        (0001-9999)
*              %D -> d    : day of gps week (0-6)
*              %H -> h    : hour code       (a=0,b=1,c=2,...,x=23)
*              %ha-> hh   : 3 hours         (00,03,06,...,21)
*              %hb-> hh   : 6 hours         (00,06,12,18)
*              %hc-> hh   : 12 hours        (00,12)
*              %t -> mm   : 15 minutes      (00,15,30,45)
*              %r -> rrrr : rover id
*              %b -> bbbb : base station id
*---------------------------------------------------------------------------------- */
int reppath(const string SrcPath,string &DstPath,gtime_t GpsTime,
	const string RovID,const string BaseID){

	gtime_t curTime;
	double eps[]={2000,1,1,0,0,0};
	int nWeek, nDow, nDoy,stat=0;
	string strRep;

	curTime.epoch2time(eps);
	DstPath=SrcPath;

	if(DstPath.find('%')==string::npos) return 0;
	if(RovID!="") stat|=repstr(DstPath,"%r",RovID);
	if(BaseID!="")stat|=repstr(DstPath,"%b",BaseID);

	if (GpsTime.time!=0){
		GpsTime.time2epoch();
		nDow=(int)floor(GpsTime.time2gpst(&nWeek)/86400.0);
		nDoy=(int)floor(GpsTime.timediff(curTime)/86400.0)+1;
		int2str(2,"0",((int)GpsTime.ep[3]/3)*3,strRep);		stat|=repstr(DstPath,"%ha",strRep);
		int2str(2,"0",((int)GpsTime.ep[3]/6)*6,strRep);		stat|=repstr(DstPath,"%hb",strRep);
		int2str(2,"0",((int)GpsTime.ep[3]/12)*12,strRep);	stat|=repstr(DstPath,"%hc",strRep);
		doul2str(4,0,"0",GpsTime.ep[0],strRep);				stat|=repstr(DstPath,"%Y",strRep);
		doul2str(2,0,"0",fmod(GpsTime.ep[0],100.0),strRep);	stat|=repstr(DstPath,"%y",strRep);
		doul2str(2,0,"0",GpsTime.ep[1],strRep);				stat|=repstr(DstPath,"%m",strRep);
		doul2str(2,0,"0",GpsTime.ep[2],strRep);				stat|=repstr(DstPath,"%d",strRep);
		doul2str(2,0,"0",GpsTime.ep[3],strRep);				stat|=repstr(DstPath,"%h",strRep);
		doul2str(2,0,"0",GpsTime.ep[4],strRep);				stat|=repstr(DstPath,"%M",strRep);
		doul2str(2,0,"0",GpsTime.ep[5],strRep);				stat|=repstr(DstPath,"%S",strRep);
		int2str(3,"0",nDoy,strRep);							stat|=repstr(DstPath,"%n",strRep);
		int2str(4,"0",nWeek,strRep);						stat|=repstr(DstPath,"%W",strRep);
		int2str(1,"0",nDow,strRep);							stat|=repstr(DstPath,"%D",strRep);
		strRep='a'+(int)GpsTime.ep[3];						stat|=repstr(DstPath,"%H",strRep);
		int2str(2,"0",((int)eps[4]/15)*15,strRep);			stat|=repstr(DstPath,"%t",strRep);
	}
	else if (DstPath.find("%ha")||DstPath.find("%hb")||DstPath.find("%hc")||DstPath.find("%Y")||
			DstPath.find("%y")||DstPath.find("%m")||DstPath.find("%d")||DstPath.find("%h")||
			DstPath.find("%M")||DstPath.find("%S")||DstPath.find("%n")||DstPath.find("%W")||
			DstPath.find("%D")||DstPath.find("%H")||DstPath.find("%t")){
		return -1; /* no valid time */
	}
	return stat;
}
/* convert double number to string ---------------------------------------------------
* 0 <= decimals <= 6 !!!
----------------------------------------------------------------------------------- */
string doul2str(int StrLen, int DecLen, const string StrFiller, const double SrcNum, string &DstStr){
	DstStr=to_string(SrcNum); /* with 6 decimal digit */
	if(DecLen>0)DstStr=DstStr.substr(0,DstStr.length()-6+DecLen); /* decimal digits */
	else DstStr=DstStr.substr(0,DstStr.length()-7);
	while (StrLen>DstStr.length()) 
		DstStr=StrFiller+DstStr;
	return DstStr;
}
/* convert int number to string --------------------------------------------------- */
string int2str(int StrLen, const string StrFiller, const int SrcNum, string &Dststr){
	Dststr=to_string(SrcNum);
	while(StrLen>Dststr.length())
		Dststr=StrFiller+Dststr;
	return Dststr;
}
/* convert string to double number ------------------------------------------------ */
int str2double(string SrcStr,double &DstNum){
	int i,fnum;
	
	if (SrcStr.length()<=0) return 0;

	for (i=0,fnum=1; i<SrcStr.length(); i++){
		if (fnum && SrcStr[i]<='9' && SrcStr[i]>='0') fnum=0;
		if (fnum==0 && SrcStr[i]=='D') { SrcStr[i]='E'; break; }
	}
	if (fnum) return 0;

	DstNum=stod(SrcStr);
	
	return 1;
}
/* convert string to int number --------------------------------------------------- */
int str2int(const string SrcStr,int &DstNum){
	int i;
	for (i=0; i<SrcStr.length(); i++){
		if (SrcStr[i]<='9'&&SrcStr[i]>='0') break;
	}
	if (i>=SrcStr.length()) return 0;

	DstNum=stoi(SrcStr);

	return 1;
}


/* time and position transfer functions ----------------------------------------------------------- */
/* get tick time ---------------------------------------------------------------------
* get current tick in ms
* args   : none
* return : current tick in ms
*---------------------------------------------------------------------------------- */
unsigned int tickget(){
#ifdef WIN32
	return (unsigned int)timeGetTime();
#else
	struct timespec tp={ 0 };
	struct timeval  tv={ 0 };

#ifdef CLOCK_MONOTONIC_RAW
	/* linux kernel > 2.6.28 */
	if (!clock_gettime(CLOCK_MONOTONIC_RAW,&tp)) {
		return tp.tv_sec*1000u+tp.tv_nsec/1000000u;
	}
	else {
		gettimeofday(&tv,NULL);
		return tv.tv_sec*1000u+tv.tv_usec/1000u;
	}
#else
	gettimeofday(&tv,NULL);
	return tv.tv_sec*1000u+tv.tv_usec/1000u;
#endif
#endif /* WIN32 */
}
/* sleep ms ----------------------------------------------------------------------- */
void sleepms(int ms){
#ifdef WIN32
	if (ms<5) Sleep(1); else Sleep(ms);
#else
	struct timespec ts;
	if (ms<=0) return;
	ts.tv_sec=(time_t)(ms/1000);
	ts.tv_nsec=(long)(ms%1000*1000000);
	nanosleep(&ts,NULL);
#endif
}
/* adjust gps week number ------------------------------------------------------------
* adjust gps week number using cpu time
* args   : int  UnWeek       I   not-adjusted gps week number
* return : adjusted gps week number
*---------------------------------------------------------------------------------- */
int adjgpsweek(int UnWeek){
	int nW;
	gtime_t curTime;
	curTime.timeget()->utc2gpst()->time2gpst(&nW);
	if (nW<1560) nW=1560; /* use 2009/12/1 if time is earlier than 2009/12/1 */
	return UnWeek+(nW-UnWeek+512)/1024*1024;
}

/* geography and astronomy functions -------------------------------------------------------------- */
/* coordinate rotation matrix ----------------------------------------------------- */
#define Rx(t,X) do { \
    (X)[0]=1.0; (X)[1]=(X)[2]=(X)[3]=(X)[6]=0.0; \
    (X)[4]=(X)[8]=cos(t); (X)[7]=sin(t); (X)[5]=-(X)[7]; \
} while (0)

#define Ry(t,X) do { \
    (X)[4]=1.0; (X)[1]=(X)[3]=(X)[5]=(X)[7]=0.0; \
    (X)[0]=(X)[8]=cos(t); (X)[2]=sin(t); (X)[6]=-(X)[2]; \
} while (0)

#define Rz(t,X) do { \
    (X)[8]=1.0; (X)[2]=(X)[5]=(X)[6]=(X)[7]=0.0; \
    (X)[0]=(X)[4]=cos(t); (X)[3]=sin(t); (X)[1]=-(X)[3]; \
} while (0)
/* astronomical arguments: f={l,l',F,D,OMG} (rad) --------------------------------- */
void ast_args(double t,double *f){
	static const double fc[][5]={ /* coefficients for iau 1980 nutation */
		{ 134.96340251, 1717915923.2178,  31.8792,  0.051635, -0.00024470 },
		{ 357.52910918,  129596581.0481,  -0.5532,  0.000136, -0.00001149 },
		{ 93.27209062, 1739527262.8478, -12.7512, -0.001037,  0.00000417 },
		{ 297.85019547, 1602961601.2090,  -6.3706,  0.006593, -0.00003169 },
		{ 125.04455501,   -6962890.2665,   7.4722,  0.007702, -0.00005939 }
	};
	double tt[4];
	int i,j;

	for (tt[0]=t,i=1; i<4; i++) tt[i]=tt[i-1]*t;
	for (i=0; i<5; i++) {
		f[i]=fc[i][0]*3600.0;
		for (j=0; j<4; j++) f[i]+=fc[i][j+1]*tt[j];
		f[i]=fmod(f[i]*AS2R,2.0*PI);
	}
}
/* sun and moon position in eci (ref [4] 5.1.1, 5.2.1) ---------------------------- */
static void sunmoonpos_eci(gtime_t tut,double *rsun,double *rmoon)
{
	const double ep2000[]={ 2000,1,1,12,0,0 };
	double t,f[5],eps,Ms,ls,rs,lm,pm,rm,sine,cose,sinp,cosp,sinl,cosl;
	gtime_t baseTime;

	baseTime.epoch2time(ep2000);
	t=tut.timediff(baseTime)/86400.0/36525.0;

	/* astronomical arguments */
	ast_args(t,f);

	/* obliquity of the ecliptic */
	eps=23.439291-0.0130042*t;
	sine=sin(eps*D2R); cose=cos(eps*D2R);

	/* sun position in eci */
	if (rsun) {
		Ms=357.5277233+35999.05034*t;
		ls=280.460+36000.770*t+1.914666471*sin(Ms*D2R)+0.019994643*sin(2.0*Ms*D2R);
		rs=AU*(1.000140612-0.016708617*cos(Ms*D2R)-0.000139589*cos(2.0*Ms*D2R));
		sinl=sin(ls*D2R); cosl=cos(ls*D2R);
		rsun[0]=rs*cosl;
		rsun[1]=rs*cose*sinl;
		rsun[2]=rs*sine*sinl;
	}
	/* moon position in eci */
	if (rmoon) {
		lm=218.32+481267.883*t+6.29*sin(f[0])-1.27*sin(f[0]-2.0*f[3])+
			0.66*sin(2.0*f[3])+0.21*sin(2.0*f[0])-0.19*sin(f[1])-0.11*sin(2.0*f[2]);
		pm=5.13*sin(f[2])+0.28*sin(f[0]+f[2])-0.28*sin(f[2]-f[0])-
			0.17*sin(f[2]-2.0*f[3]);
		rm=RE_WGS84/sin((0.9508+0.0518*cos(f[0])+0.0095*cos(f[0]-2.0*f[3])+
			0.0078*cos(2.0*f[3])+0.0028*cos(2.0*f[0]))*D2R);
		sinl=sin(lm*D2R); cosl=cos(lm*D2R);
		sinp=sin(pm*D2R); cosp=cos(pm*D2R);
		rmoon[0]=rm*cosp*cosl;
		rmoon[1]=rm*(cose*cosp*sinl-sine*sinp);
		rmoon[2]=rm*(sine*cosp*sinl+cose*sinp);
	}
}
/* iau 1980 nutation -------------------------------------------------------------- */
static void nut_iau1980(double t,const double *f,double *dpsi,double *deps)
{
	static const double nut[106][10]={
		{ 0,   0,   0,   0,   1, -6798.4, -171996, -174.2, 92025,   8.9 },
		{ 0,   0,   2,  -2,   2,   182.6,  -13187,   -1.6,  5736,  -3.1 },
		{ 0,   0,   2,   0,   2,    13.7,   -2274,   -0.2,   977,  -0.5 },
		{ 0,   0,   0,   0,   2, -3399.2,    2062,    0.2,  -895,   0.5 },
		{ 0,  -1,   0,   0,   0,  -365.3,   -1426,    3.4,    54,  -0.1 },
		{ 1,   0,   0,   0,   0,    27.6,     712,    0.1,    -7,   0.0 },
		{ 0,   1,   2,  -2,   2,   121.7,    -517,    1.2,   224,  -0.6 },
		{ 0,   0,   2,   0,   1,    13.6,    -386,   -0.4,   200,   0.0 },
		{ 1,   0,   2,   0,   2,     9.1,    -301,    0.0,   129,  -0.1 },
		{ 0,  -1,   2,  -2,   2,   365.2,     217,   -0.5,   -95,   0.3 },
		{ -1,   0,   0,   2,   0,    31.8,     158,    0.0,    -1,   0.0 },
		{ 0,   0,   2,  -2,   1,   177.8,     129,    0.1,   -70,   0.0 },
		{ -1,   0,   2,   0,   2,    27.1,     123,    0.0,   -53,   0.0 },
		{ 1,   0,   0,   0,   1,    27.7,      63,    0.1,   -33,   0.0 },
		{ 0,   0,   0,   2,   0,    14.8,      63,    0.0,    -2,   0.0 },
		{ -1,   0,   2,   2,   2,     9.6,     -59,    0.0,    26,   0.0 },
		{ -1,   0,   0,   0,   1,   -27.4,     -58,   -0.1,    32,   0.0 },
		{ 1,   0,   2,   0,   1,     9.1,     -51,    0.0,    27,   0.0 },
		{ -2,   0,   0,   2,   0,  -205.9,     -48,    0.0,     1,   0.0 },
		{ -2,   0,   2,   0,   1,  1305.5,      46,    0.0,   -24,   0.0 },
		{ 0,   0,   2,   2,   2,     7.1,     -38,    0.0,    16,   0.0 },
		{ 2,   0,   2,   0,   2,     6.9,     -31,    0.0,    13,   0.0 },
		{ 2,   0,   0,   0,   0,    13.8,      29,    0.0,    -1,   0.0 },
		{ 1,   0,   2,  -2,   2,    23.9,      29,    0.0,   -12,   0.0 },
		{ 0,   0,   2,   0,   0,    13.6,      26,    0.0,    -1,   0.0 },
		{ 0,   0,   2,  -2,   0,   173.3,     -22,    0.0,     0,   0.0 },
		{ -1,   0,   2,   0,   1,    27.0,      21,    0.0,   -10,   0.0 },
		{ 0,   2,   0,   0,   0,   182.6,      17,   -0.1,     0,   0.0 },
		{ 0,   2,   2,  -2,   2,    91.3,     -16,    0.1,     7,   0.0 },
		{ -1,   0,   0,   2,   1,    32.0,      16,    0.0,    -8,   0.0 },
		{ 0,   1,   0,   0,   1,   386.0,     -15,    0.0,     9,   0.0 },
		{ 1,   0,   0,  -2,   1,   -31.7,     -13,    0.0,     7,   0.0 },
		{ 0,  -1,   0,   0,   1,  -346.6,     -12,    0.0,     6,   0.0 },
		{ 2,   0,  -2,   0,   0, -1095.2,      11,    0.0,     0,   0.0 },
		{ -1,   0,   2,   2,   1,     9.5,     -10,    0.0,     5,   0.0 },
		{ 1,   0,   2,   2,   2,     5.6,      -8,    0.0,     3,   0.0 },
		{ 0,  -1,   2,   0,   2,    14.2,      -7,    0.0,     3,   0.0 },
		{ 0,   0,   2,   2,   1,     7.1,      -7,    0.0,     3,   0.0 },
		{ 1,   1,   0,  -2,   0,   -34.8,      -7,    0.0,     0,   0.0 },
		{ 0,   1,   2,   0,   2,    13.2,       7,    0.0,    -3,   0.0 },
		{ -2,   0,   0,   2,   1,  -199.8,      -6,    0.0,     3,   0.0 },
		{ 0,   0,   0,   2,   1,    14.8,      -6,    0.0,     3,   0.0 },
		{ 2,   0,   2,  -2,   2,    12.8,       6,    0.0,    -3,   0.0 },
		{ 1,   0,   0,   2,   0,     9.6,       6,    0.0,     0,   0.0 },
		{ 1,   0,   2,  -2,   1,    23.9,       6,    0.0,    -3,   0.0 },
		{ 0,   0,   0,  -2,   1,   -14.7,      -5,    0.0,     3,   0.0 },
		{ 0,  -1,   2,  -2,   1,   346.6,      -5,    0.0,     3,   0.0 },
		{ 2,   0,   2,   0,   1,     6.9,      -5,    0.0,     3,   0.0 },
		{ 1,  -1,   0,   0,   0,    29.8,       5,    0.0,     0,   0.0 },
		{ 1,   0,   0,  -1,   0,   411.8,      -4,    0.0,     0,   0.0 },
		{ 0,   0,   0,   1,   0,    29.5,      -4,    0.0,     0,   0.0 },
		{ 0,   1,   0,  -2,   0,   -15.4,      -4,    0.0,     0,   0.0 },
		{ 1,   0,  -2,   0,   0,   -26.9,       4,    0.0,     0,   0.0 },
		{ 2,   0,   0,  -2,   1,   212.3,       4,    0.0,    -2,   0.0 },
		{ 0,   1,   2,  -2,   1,   119.6,       4,    0.0,    -2,   0.0 },
		{ 1,   1,   0,   0,   0,    25.6,      -3,    0.0,     0,   0.0 },
		{ 1,  -1,   0,  -1,   0, -3232.9,      -3,    0.0,     0,   0.0 },
		{ -1,  -1,   2,   2,   2,     9.8,      -3,    0.0,     1,   0.0 },
		{ 0,  -1,   2,   2,   2,     7.2,      -3,    0.0,     1,   0.0 },
		{ 1,  -1,   2,   0,   2,     9.4,      -3,    0.0,     1,   0.0 },
		{ 3,   0,   2,   0,   2,     5.5,      -3,    0.0,     1,   0.0 },
		{ -2,   0,   2,   0,   2,  1615.7,      -3,    0.0,     1,   0.0 },
		{ 1,   0,   2,   0,   0,     9.1,       3,    0.0,     0,   0.0 },
		{ -1,   0,   2,   4,   2,     5.8,      -2,    0.0,     1,   0.0 },
		{ 1,   0,   0,   0,   2,    27.8,      -2,    0.0,     1,   0.0 },
		{ -1,   0,   2,  -2,   1,   -32.6,      -2,    0.0,     1,   0.0 },
		{ 0,  -2,   2,  -2,   1,  6786.3,      -2,    0.0,     1,   0.0 },
		{ -2,   0,   0,   0,   1,   -13.7,      -2,    0.0,     1,   0.0 },
		{ 2,   0,   0,   0,   1,    13.8,       2,    0.0,    -1,   0.0 },
		{ 3,   0,   0,   0,   0,     9.2,       2,    0.0,     0,   0.0 },
		{ 1,   1,   2,   0,   2,     8.9,       2,    0.0,    -1,   0.0 },
		{ 0,   0,   2,   1,   2,     9.3,       2,    0.0,    -1,   0.0 },
		{ 1,   0,   0,   2,   1,     9.6,      -1,    0.0,     0,   0.0 },
		{ 1,   0,   2,   2,   1,     5.6,      -1,    0.0,     1,   0.0 },
		{ 1,   1,   0,  -2,   1,   -34.7,      -1,    0.0,     0,   0.0 },
		{ 0,   1,   0,   2,   0,    14.2,      -1,    0.0,     0,   0.0 },
		{ 0,   1,   2,  -2,   0,   117.5,      -1,    0.0,     0,   0.0 },
		{ 0,   1,  -2,   2,   0,  -329.8,      -1,    0.0,     0,   0.0 },
		{ 1,   0,  -2,   2,   0,    23.8,      -1,    0.0,     0,   0.0 },
		{ 1,   0,  -2,  -2,   0,    -9.5,      -1,    0.0,     0,   0.0 },
		{ 1,   0,   2,  -2,   0,    32.8,      -1,    0.0,     0,   0.0 },
		{ 1,   0,   0,  -4,   0,   -10.1,      -1,    0.0,     0,   0.0 },
		{ 2,   0,   0,  -4,   0,   -15.9,      -1,    0.0,     0,   0.0 },
		{ 0,   0,   2,   4,   2,     4.8,      -1,    0.0,     0,   0.0 },
		{ 0,   0,   2,  -1,   2,    25.4,      -1,    0.0,     0,   0.0 },
		{ -2,   0,   2,   4,   2,     7.3,      -1,    0.0,     1,   0.0 },
		{ 2,   0,   2,   2,   2,     4.7,      -1,    0.0,     0,   0.0 },
		{ 0,  -1,   2,   0,   1,    14.2,      -1,    0.0,     0,   0.0 },
		{ 0,   0,  -2,   0,   1,   -13.6,      -1,    0.0,     0,   0.0 },
		{ 0,   0,   4,  -2,   2,    12.7,       1,    0.0,     0,   0.0 },
		{ 0,   1,   0,   0,   2,   409.2,       1,    0.0,     0,   0.0 },
		{ 1,   1,   2,  -2,   2,    22.5,       1,    0.0,    -1,   0.0 },
		{ 3,   0,   2,  -2,   2,     8.7,       1,    0.0,     0,   0.0 },
		{ -2,   0,   2,   2,   2,    14.6,       1,    0.0,    -1,   0.0 },
		{ -1,   0,   0,   0,   2,   -27.3,       1,    0.0,    -1,   0.0 },
		{ 0,   0,  -2,   2,   1,  -169.0,       1,    0.0,     0,   0.0 },
		{ 0,   1,   2,   0,   1,    13.1,       1,    0.0,     0,   0.0 },
		{ -1,   0,   4,   0,   2,     9.1,       1,    0.0,     0,   0.0 },
		{ 2,   1,   0,  -2,   0,   131.7,       1,    0.0,     0,   0.0 },
		{ 2,   0,   0,   2,   0,     7.1,       1,    0.0,     0,   0.0 },
		{ 2,   0,   2,  -2,   1,    12.8,       1,    0.0,    -1,   0.0 },
		{ 2,   0,  -2,   0,   1,  -943.2,       1,    0.0,     0,   0.0 },
		{ 1,  -1,   0,  -2,   0,   -29.3,       1,    0.0,     0,   0.0 },
		{ -1,   0,   0,   1,   1,  -388.3,       1,    0.0,     0,   0.0 },
		{ -1,  -1,   0,   2,   1,    35.0,       1,    0.0,     0,   0.0 },
		{ 0,   1,   0,   1,   0,    27.3,       1,    0.0,     0,   0.0 }
	};
	double ang;
	int i,j;

	*dpsi=*deps=0.0;

	for (i=0; i<106; i++) {
		ang=0.0;
		for (j=0; j<5; j++) ang+=nut[i][j]*f[j];
		*dpsi+=(nut[i][6]+nut[i][7]*t)*sin(ang);
		*deps+=(nut[i][8]+nut[i][9]*t)*cos(ang);
	}
	*dpsi*=1E-4*AS2R; /* 0.1 mas -> rad */
	*deps*=1E-4*AS2R;
}
/* eci to ecef transformation matrix -------------------------------------------
* compute eci to ecef transformation matrix
* args   : gtime_t tutc     I   time in utc
*          double *erpv     I   erp values {xp,yp,ut1_utc,lod} (rad,rad,s,s/d)
*          double *U        O   eci to ecef transformation matrix (3 x 3)
*          double *gmst     IO  greenwich mean sidereal time (rad)
*                               (NULL: no output)
* return : none
* note   : see ref [3] chap 5
*          not thread-safe
*-----------------------------------------------------------------------------*/
void eci2ecef(gtime_t tutc,const double *erpv,double *U,double *gmst)
{
	const double ep2000[]={ 2000,1,1,12,0,0 };
	static gtime_t tutc_;
	static double U_[9],gmst_;
	gtime_t tgps,baseTime;
	double eps,ze,th,z,t,t2,t3,dpsi,deps,gast,f[5];
	double R1[9],R2[9],R3[9],R[9],W[9],N[9],P[9],NP[9];
	int i;

	if (fabs(tutc.timediff(tutc_))<0.01) { /* read cache */
		for (i=0; i<9; i++) U[i]=U_[i];
		if (gmst) *gmst=gmst_;
		return;
	}
	tutc_=tutc;

	/* terrestrial time */
	tgps=tutc_;
	tgps.utc2gpst();
	baseTime.epoch2time(ep2000);
	t=(tgps.timediff(baseTime)+19.0+32.184)/86400.0/36525.0;
	t2=t*t; t3=t2*t;

	/* astronomical arguments */
	ast_args(t,f);

	/* iau 1976 precession */
	ze=(2306.2181*t+0.30188*t2+0.017998*t3)*AS2R;
	th=(2004.3109*t-0.42665*t2-0.041833*t3)*AS2R;
	z =(2306.2181*t+1.09468*t2+0.018203*t3)*AS2R;
	eps=(84381.448-46.8150*t-0.00059*t2+0.001813*t3)*AS2R;
	Rz(-z,R1); Ry(th,R2); Rz(-ze,R3);
	matmul_pnt("NN",3,3,3,1.0,R1,R2,0.0,R);
	matmul_pnt("NN",3,3,3,1.0,R,R3,0.0,P); /* P=Rz(-z)*Ry(th)*Rz(-ze) */

									   /* iau 1980 nutation */
	nut_iau1980(t,f,&dpsi,&deps);
	Rx(-eps-deps,R1); Rz(-dpsi,R2); Rx(eps,R3);
	matmul_pnt("NN",3,3,3,1.0,R1,R2,0.0,R);
	matmul_pnt("NN",3,3,3,1.0,R,R3,0.0,N); /* N=Rx(-eps)*Rz(-dspi)*Rx(eps) */

									   /* greenwich aparent sidereal time (rad) */
	gmst_=tutc_.utc2gmst(erpv[2]);
	gast=gmst_+dpsi*cos(eps);
	gast+=(0.00264*sin(f[4])+0.000063*sin(2.0*f[4]))*AS2R;

	/* eci to ecef transformation matrix */
	Ry(-erpv[0],R1); Rx(-erpv[1],R2); Rz(gast,R3);
	matmul_pnt("NN",3,3,3,1.0,R1,R2,0.0,W);
	matmul_pnt("NN",3,3,3,1.0,W,R3,0.0,R); /* W=Ry(-xp)*Rx(-yp) */
	matmul_pnt("NN",3,3,3,1.0,N,P,0.0,NP);
	matmul_pnt("NN",3,3,3,1.0,R,NP,0.0,U_); /* U=W*Rz(gast)*N*P */

	for (i=0; i<9; i++) U[i]=U_[i];
	if (gmst) *gmst=gmst_;
}
/* sun and moon position -------------------------------------------------------
* get sun and moon position in ecef
* args   : gtime_t tut      I   time in ut1
*          double *erpv     I   erp value {xp,yp,ut1_utc,lod} (rad,rad,s,s/d)
*          double *rsun     IO  sun position in ecef  (m) (NULL: not output)
*          double *rmoon    IO  moon position in ecef (m) (NULL: not output)
*          double *gmst     O   gmst (rad)
* return : none
*-----------------------------------------------------------------------------*/
void sunmoonpos(gtime_t UT1Time,double *ERPValue,double *SunPos,
	double *MoonPos,double *gmst)
{
	double rs[3],rm[3],U[9],gmst_;

	UT1Time.timeadd(ERPValue[2]); /* utc -> ut1 */

						   /* sun and moon position in eci */
	sunmoonpos_eci(UT1Time,SunPos ? rs : NULL,MoonPos ? rm : NULL);

	/* eci to ecef transformation matrix */
	eci2ecef(UT1Time,ERPValue,U,&gmst_);

	/* sun and moon postion in ecef */
	if (SunPos) matmul_pnt("NN",3,1,3,1.0,U,rs,0.0,SunPos);
	if (MoonPos) matmul_pnt("NN",3,1,3,1.0,U,rm,0.0,MoonPos);
	if (gmst) *gmst=gmst_;
}
/* get earth rotation parameter values -------------------------------------------- */
int geterpv(const erp_t *Earth_Par,gtime_t GPS_Time,double *Earth_Value){
	const double ep[]={ 2000,1,1,12,0,0 };
	const gtime_t time(ep);
	double mjd,day,a;
	int i,j,k;

	if (Earth_Par->n<=0) return 0;

	mjd=51544.5+GPS_Time.gpst2utc()->timediff(time)/86400.0;

	if (mjd<=Earth_Par->data[0].mjd) {
		day=mjd-Earth_Par->data[0].mjd;
		Earth_Value[0]=Earth_Par->data[0].xp     +Earth_Par->data[0].xpr*day;
		Earth_Value[1]=Earth_Par->data[0].yp     +Earth_Par->data[0].ypr*day;
		Earth_Value[2]=Earth_Par->data[0].ut1_utc-Earth_Par->data[0].lod*day;
		Earth_Value[3]=Earth_Par->data[0].lod;
		return 1;
	}
	if (mjd>=Earth_Par->data[Earth_Par->n-1].mjd) {
		day=mjd-Earth_Par->data[Earth_Par->n-1].mjd;
		Earth_Value[0]=Earth_Par->data[Earth_Par->n-1].xp     +Earth_Par->data[Earth_Par->n-1].xpr*day;
		Earth_Value[1]=Earth_Par->data[Earth_Par->n-1].yp     +Earth_Par->data[Earth_Par->n-1].ypr*day;
		Earth_Value[2]=Earth_Par->data[Earth_Par->n-1].ut1_utc-Earth_Par->data[Earth_Par->n-1].lod*day;
		Earth_Value[3]=Earth_Par->data[Earth_Par->n-1].lod;
		return 1;
	}
	for (j=0,k=Earth_Par->n-1; j<k-1;) {
		i=(j+k)/2;
		if (mjd<Earth_Par->data[i].mjd) k=i; else j=i;
	}
	if (Earth_Par->data[j].mjd==Earth_Par->data[j+1].mjd) {
		a=0.5;
	}
	else {
		a=(mjd-Earth_Par->data[j].mjd)/(Earth_Par->data[j+1].mjd-Earth_Par->data[j].mjd);
	}
	Earth_Value[0]=(1.0-a)*Earth_Par->data[j].xp     +a*Earth_Par->data[j+1].xp;
	Earth_Value[1]=(1.0-a)*Earth_Par->data[j].yp     +a*Earth_Par->data[j+1].yp;
	Earth_Value[2]=(1.0-a)*Earth_Par->data[j].ut1_utc+a*Earth_Par->data[j+1].ut1_utc;
	Earth_Value[3]=(1.0-a)*Earth_Par->data[j].lod    +a*Earth_Par->data[j+1].lod;
	return 1;
}

/* satellite data functions ----------------------------------------------------------------------- */
/* satellite system+prn/slot number to satellite number ------------------------------
* convert satellite system+prn/slot number to satellite number
* args   : int    SatSys       I   satellite system (SYS_GPS,SYS_GLO,...)
*          int    PrnNum       I   satellite prn/slot number
* return : satellite number (0:error)
*---------------------------------------------------------------------------------- */
int satno(int SatSys, int PrnNum)
{
    if (PrnNum<=0) return 0;
    switch (SatSys) {
        case SYS_GPS:
            if (PrnNum<MINPRNGPS||MAXPRNGPS<PrnNum) return 0;
            return PrnNum-MINPRNGPS+1;
        case SYS_GLO:
            if (PrnNum<MINPRNGLO||MAXPRNGLO<PrnNum) return 0;
            return NSATGPS+PrnNum-MINPRNGLO+1;
        case SYS_GAL:
            if (PrnNum<MINPRNGAL||MAXPRNGAL<PrnNum) return 0;
            return NSATGPS+NSATGLO+PrnNum-MINPRNGAL+1;
        case SYS_QZS:
            if (PrnNum<MINPRNQZS||MAXPRNQZS<PrnNum) return 0;
            return NSATGPS+NSATGLO+NSATGAL+PrnNum-MINPRNQZS+1;
        case SYS_CMP:
            if (PrnNum<MINPRNCMP||MAXPRNCMP<PrnNum) return 0;
            return NSATGPS+NSATGLO+NSATGAL+NSATQZS+PrnNum-MINPRNCMP+1;
        case SYS_IRN:
            if (PrnNum<MINPRNIRN||MAXPRNIRN<PrnNum) return 0;
            return NSATGPS+NSATGLO+NSATGAL+NSATQZS+NSATCMP+PrnNum-MINPRNIRN+1;
        case SYS_LEO:
            if (PrnNum<MINPRNLEO||MAXPRNLEO<PrnNum) return 0;
            return NSATGPS+NSATGLO+NSATGAL+NSATQZS+NSATCMP+NSATIRN+
                   PrnNum-MINPRNLEO+1;
        case SYS_SBS:
            if (PrnNum<MINPRNSBS||MAXPRNSBS<PrnNum) return 0;
            return NSATGPS+NSATGLO+NSATGAL+NSATQZS+NSATCMP+NSATIRN+NSATLEO+
                   PrnNum-MINPRNSBS+1;
    }
    return 0;
}
/* satellite number to satellite system ----------------------------------------------
* convert satellite number to satellite system
* args   : int    SatNum       I   satellite number (1-MAXSAT)
*          int    *PrnNum      IO  satellite prn/slot number (NULL: no output)
* return : satellite system (SYS_GPS,SYS_GLO,...)
*---------------------------------------------------------------------------------- */
int satsys(int SatNum, int *PrnNum)
{
    int satSys=SYS_NONE;
    if (SatNum<=0||MAXSAT<SatNum) SatNum=0;
    else if (SatNum<=NSATGPS) {
        satSys=SYS_GPS; SatNum+=MINPRNGPS-1;
    }
    else if ((SatNum-=NSATGPS)<=NSATGLO) {
        satSys=SYS_GLO; SatNum+=MINPRNGLO-1;
    }
    else if ((SatNum-=NSATGLO)<=NSATGAL) {
        satSys=SYS_GAL; SatNum+=MINPRNGAL-1;
    }
    else if ((SatNum-=NSATGAL)<=NSATQZS) {
        satSys=SYS_QZS; SatNum+=MINPRNQZS-1; 
    }
    else if ((SatNum-=NSATQZS)<=NSATCMP) {
        satSys=SYS_CMP; SatNum+=MINPRNCMP-1; 
    }
    else if ((SatNum-=NSATCMP)<=NSATIRN) {
		satSys=SYS_IRN; SatNum+=MINPRNIRN-1;
    }
    else if ((SatNum-=NSATIRN)<=NSATLEO) {
		satSys=SYS_LEO; SatNum+=MINPRNLEO-1;
    }
    else if ((SatNum-=NSATLEO)<=NSATSBS) {
		satSys=SYS_SBS; SatNum+=MINPRNSBS-1;
    }
    else SatNum=0;
    if (PrnNum) *PrnNum=SatNum;
    return satSys;
}
/* satellite id to satellite number --------------------------------------------------
* convert satellite id to satellite number
* args   : string   SatID       I   satellite id (nn,Gnn,Rnn,Enn,Jnn,Cnn,Inn or Snn)
* return : satellite number     (0: error)
* notes  : 120-142 and 193-199 are also recognized as sbas and qzss
*---------------------------------------------------------------------------------- */
int satid2no(string SatID){
	int satSys,prnNum;
	char code;
	string aaa=" 0123456789";
	
	if (aaa.find(SatID[0])!=string::npos) {
		if (str2int((SatID.substr(1,2)),prnNum)==0) return 0;
		if      (MINPRNGPS<=prnNum&&prnNum<=MAXPRNGPS) satSys=SYS_GPS;
		else if (MINPRNSBS<=prnNum&&prnNum<=MAXPRNSBS) satSys=SYS_SBS;
		else if (MINPRNQZS<=prnNum&&prnNum<=MAXPRNQZS) satSys=SYS_QZS;
		else return 0;
		return satno(satSys,prnNum);
	}
	
	if(str2int((SatID.substr(1,2)),prnNum)==0) return 0;
	code=SatID[0];

	switch (code) {
		case 'G': satSys=SYS_GPS; prnNum+=MINPRNGPS-1; break;
		case 'R': satSys=SYS_GLO; prnNum+=MINPRNGLO-1; break;
		case 'E': satSys=SYS_GAL; prnNum+=MINPRNGAL-1; break;
		case 'J': satSys=SYS_QZS; prnNum+=MINPRNQZS-1; break;
		case 'C': satSys=SYS_CMP; prnNum+=MINPRNCMP-1; break;
		case 'I': satSys=SYS_IRN; prnNum+=MINPRNIRN-1; break;
		case 'L': satSys=SYS_LEO; prnNum+=MINPRNLEO-1; break;
		case 'S': satSys=SYS_SBS; prnNum+=MINPRNSBS-1; break;
		default: return 0;
	}
	return satno(satSys,prnNum);
}
/* satellite number to satellite id ----------------------------------------------- */
int satno2id(int sat,string &SatID){
	int prn;
	string buf;
	int sys=satsys(sat,&prn);
	switch (sys){
		case SYS_GPS: SatID="G"+int2str(2,"0",prn-MINPRNGPS+1,buf); return SYS_GPS;
		case SYS_GLO: SatID="R"+int2str(2,"0",prn-MINPRNGLO+1,buf); return SYS_GLO;
		case SYS_GAL: SatID="E"+int2str(2,"0",prn-MINPRNGAL+1,buf); return SYS_GAL;
		case SYS_QZS: SatID="J"+int2str(2,"0",prn-MINPRNQZS+1,buf); return SYS_QZS;
		case SYS_CMP: SatID="C"+int2str(2,"0",prn-MINPRNCMP+1,buf); return SYS_CMP;
		case SYS_IRN: SatID="I"+int2str(2,"0",prn-MINPRNIRN+1,buf); return SYS_IRN;
		case SYS_LEO: SatID="L"+int2str(2,"0",prn-MINPRNLEO+1,buf); return SYS_LEO;
		case SYS_SBS: SatID="S"+int2str(2,"0",prn-MINPRNSBS+1,buf); return SYS_SBS;
	}
	SatID="";
	return SYS_NONE;
}
/* satellite carrier wave length -------------------------------------------------- */
double satwavelen(int SatNum,int FrqNum,const nav_t *NavData){
	const double freq_glos[]={ FREQ1_GLO,FREQ2_GLO };
	const double dfrq_glos[]={ DFRQ1_GLO,DFRQ2_GLO };
	int prn;
	int i,satSys=satsys(SatNum,&prn);

	if (satSys==SYS_GLO) {
		if (0<=FrqNum&&FrqNum<=1) {
			for (i=0; i<NavData->ng; i++) {
				if (NavData->geph[i].sat!=SatNum) continue;
				return CLIGHT/(freq_glos[FrqNum]+dfrq_glos[FrqNum]*NavData->geph[i].frq);
			}
		}
		else if (FrqNum==2) { /* L3 */
			return CLIGHT/FREQ3_GLO;
		}
	}
	else if (satSys==SYS_CMP) {
		/* BeiDou-2 satellite */
		if (prn<18) {
			if (FrqNum==0) return CLIGHT/FREQ1_CMP; /* B1 */
			else if (FrqNum==1) return CLIGHT/FREQ2_CMP; /* B2 */
			else if (FrqNum==2) return CLIGHT/FREQ3_CMP; /* B3 */
		}
		/* BeiDou-3 satellite */
		else {
			if (FrqNum==0) return CLIGHT/FREQ1_CMP_3; /* B1C */
			else if (FrqNum==1) return CLIGHT/FREQ2_CMP_3; /* B2a */
			else if (FrqNum==2) return CLIGHT/FREQ3_CMP; /* B3 */
		}
	}
	else {
		if (FrqNum==0) return CLIGHT/FREQ1; /* L1/E1 */
		else if (FrqNum==1) return CLIGHT/FREQ2; /* L2 */
		else if (FrqNum==2) return CLIGHT/FREQ5; /* L5/E5a */
		else if (FrqNum==3) return CLIGHT/FREQ6; /* L6/LEX */
		else if (FrqNum==4) return CLIGHT/FREQ7; /* E5b */
		else if (FrqNum==5) return CLIGHT/FREQ8; /* E5a+b */
		else if (FrqNum==6) return CLIGHT/FREQ9; /* S */
	}
	return 0.0;
}

/* observation and code transfer functions -------------------------------------------------------- */
/* obs type string to obs code -------------------------------------------------------
* convert obs code type string to obs code
* args   : string   ObsCode   I     obs code string ("1C","1P","1Y",...)
*          int      *ObsFre  IO     frequency (1:L1,2:L2,3:L5,4:L6,5:L7,6:L8,0:err)
*                               (NULL: no output)
* return : obs code (CODE_???)
* notes  : obs codes are based on reference [6] and qzss extension
*---------------------------------------------------------------------------------- */
unsigned char obs2code(string ObsCode, int *ObsFre)
{
    int i;
    if (ObsFre) *ObsFre=0;
    for (i=1;i<60;i++) {
        if (ObsCodes[i].compare(ObsCode)!=0) continue;
        if (ObsFre) *ObsFre=ObsFreqs[i];
        return (unsigned char)i;
    }
    return CODE_NONE;
}
/* obs code to obs code string -------------------------------------------------------
* convert obs code to obs code string
* args   : unsigned char code I obs code (CODE_???)
*          int    *freq  IO     frequency (NULL: no output)
*                               (1:L1/E1, 2:L2/B1, 3:L5/E5a/L3, 4:L6/LEX/B3,
                                 5:E5b/B2, 6:E5(a+b), 7:S)
* return : obs code string ("1C","1P","1P",...)
* notes  : obs codes are based on reference [6] and qzss extension
*---------------------------------------------------------------------------------- */
string code2obs(unsigned char ObsCode, int *ObsFre)
{
    if (ObsFre) *ObsFre=0;
    if (ObsCode<=CODE_NONE||MAXCODE<ObsCode) return "";
    if (ObsFre) *ObsFre=ObsFreqs[ObsCode];
    return ObsCodes[ObsCode];
}
/* satellite code to satellite system --------------------------------------------- */
int code2sys(char SysCode) {
	if (SysCode=='G'||SysCode==' ') return SYS_GPS;
    if (SysCode=='R') return SYS_GLO;
    if (SysCode=='E') return SYS_GAL; /* extension to sp3-c */
    if (SysCode=='J') return SYS_QZS; /* extension to sp3-c */
    if (SysCode=='C') return SYS_CMP; /* extension to sp3-c */
    if (SysCode=='L') return SYS_LEO; /* extension to sp3-c */
    return SYS_NONE;
}
/* get code priority -----------------------------------------------------------------
* get code priority for multiple codes in a frequency
* args   : int    sys     I     system (SYS_???)
*          unsigned char code I obs code (CODE_???)
*          string opt    I      code options (NULL:no option)
* return : priority (15:highest-1:lowest,0:error)
*---------------------------------------------------------------------------------- */
int getcodepri(int SatSys, unsigned char ObsCode, string CodeOpt)
{
	if (ObsCode==0) return 0;
    size_t strPos;
    string strOpt;
    string strObs,strBuff;
    int i,j;
    
    switch (SatSys) {
        case SYS_GPS: i=0; strOpt="-GL"; break;
        case SYS_GLO: i=1; strOpt="-RL"; break;
        case SYS_GAL: i=2; strOpt="-EL"; break;
        case SYS_QZS: i=3; strOpt="-JL"; break;
        case SYS_SBS: i=4; strOpt="-SL"; break;
        case SYS_CMP: i=5; strOpt="-CL"; break;
        case SYS_IRN: i=6; strOpt="-IL"; break;
        default: return 0;
    }
    strObs=code2obs(ObsCode,&j);
    
    /* parse code options */
    for (strPos=0;strPos<CodeOpt.length();strPos++) {
        if (CodeOpt.substr(strPos,3)==strOpt){
        	strBuff=CodeOpt.substr(strPos+3,2);
        	if (strBuff[0]!=strObs[0]) continue;
        	return strBuff[1]==strObs[1]?15:0;
    	}
    }
    /* search code priority */
    return ((strPos=CodePriors[i][j-1].find(strObs[1]))!=string::npos)?14-(int)strPos:0;
}

/* system functions ------------------------------------------------------------------------------- */
/* add fatal callback function -------------------------------------------------------
* add fatal callback function for mat(),zeros(),imat()
* args   : fatalfunc_t *func I  callback function
* return : none
* notes  : if malloc() failed in return : none
*---------------------------------------------------------------------------------- */
void add_fatal(fatalfunc_t *func)
{
	FatalFunc = func;
}

/* GPS data functions ----------------------------------------------------------------------------- */
/* compare observation data ------------------------------------------------------- */
int cmpobs(const void *SrcObs1,const void *SrcObs2){
	obsd_t *pObs1=(obsd_t *)SrcObs1,*pObs2=(obsd_t *)SrcObs2;
	double diffTime=pObs1->time.timediff(pObs2->time);
	if (fabs(diffTime)>DTTOL) return diffTime<0?-1:1;
	if (pObs1->rcv!=pObs2->rcv) return (int)pObs1->rcv-(int)pObs2->rcv;
	return (int)pObs1->sat-(int)pObs2->sat;
}
/* arrange observation data ------------------------------------------------------- */
int sortobs(obs_t &SrcObs){
	int i,j,n;

	if (SrcObs.n<=0) return 0;

	qsort(&SrcObs.data[0],SrcObs.n,sizeof(obsd_t),cmpobs);

	/* delete duplicated data */
	for (i=j=0; i<SrcObs.n; i++) {
		if (SrcObs.data[i].sat!=SrcObs.data[j].sat||
			SrcObs.data[i].rcv!=SrcObs.data[j].rcv||
			SrcObs.data[i].time.timediff(SrcObs.data[j].time)!=0.0) {
			SrcObs.data[++j]=SrcObs.data[i];
		}
	}
	SrcObs.n=j+1;

	for (i=n=0; i<SrcObs.n; i=j,n++) {
		for (j=i+1; j<SrcObs.n; j++) {
			if (SrcObs.data[j].time.timediff(SrcObs.data[i].time)>DTTOL) break;
		}
	}
	return n;
}
/* station-cross doppler single-difference ---------------------------------------- */
double dopsingle_d(const obsd_t *rov,const obsd_t *bas,const int fre){
	if (bas){
		double drov=rov->D[fre];
		double dbas=bas->D[fre];
		return drov==0.0||dbas==0.0 ? 0.0 : drov-dbas;
	}
	else return rov->D[fre];
}
/* station-cross single-difference observation ---------------------------------------
* argv   : int  fre frequency type (<NFREQ:L,>NFREQ:P)
* --------------------------------------------------------------------------------- */
double single_diff(const obsd_t *rov,const obsd_t *bas,const int fre){
	if (bas){
		double prov=fre<NFREQ ? rov->L[fre] : rov->P[fre-NFREQ];
		double pbas=fre<NFREQ ? bas->L[fre] : bas->P[fre-NFREQ];
		return prov==0.0||pbas==0.0 ? 0.0 : prov-pbas;
	}
	else 
		return fre<NFREQ ? rov->L[fre] : rov->P[fre-NFREQ];
}
/* compute geometry-free combination -------------------------------------------------
* argv   : int  fres    frequency combination 
*                       (1:L1L2,2:L1L5;3:L2L5,4:P1P2,5:P1P5,6:P2P5)
* if *bas==NULL, return rov geometry-free combination
* --------------------------------------------------------------------------------- */
double geometry_free(const int fres,const obsd_t *rov,const obsd_t *bas,const double *lam){
	int f1=0,f2=0,fre=0;
	double p1=0.0,p2=0.0;

	if (fres>NFREQ) fre=NFREQ;
	switch (fres%3){
		case 1: f1=0; f2=1; break;
		case 2: f1=0; f2=2; break;
		case 0: f1=1; f2=2; break;
		default: f1=0; f2=1;
	}

	p1=single_diff(rov,bas,f1+fre) * (fre==NFREQ ? 1.0 : lam[f1]);
	p2=single_diff(rov,bas,f2+fre) * (fre==NFREQ ? 1.0 : lam[f2]);
	
	return p1==0.0||p2==0.0?0.0:p1-p2;
}
/* compute ionosphere-free combination -----------------------------------------------
* argv   : int  fres    frequency combination 
*                       (1:L1L2,2:L1L5;3:L2L5,4:P1P2,5:P1P5,6:P2P5)
* --------------------------------------------------------------------------------- */
double iono_free(const int fres,const obsd_t *rov,const obsd_t *bas,const double *lam){
	int f1=0,f2=0,fre=0;
	double p1=0.0,p2=0.0;

	if (fres>NFREQ) fre=NFREQ;
	switch (fres%3){
		case 1: f1=0; f2=1; break;
		case 2: f1=0; f2=2; break;
		case 0: f1=1; f2=2; break;
		default: f1=0; f2=1;
	}
	double gamma=SQR(lam[f2])/SQR(lam[f1]); /* f1^2/f2^2 */

	p1=single_diff(rov,bas,f1+fre) * (fre==NFREQ ? 1.0 : lam[f1]);
	p2=single_diff(rov,bas,f2+fre) * (fre==NFREQ ? 1.0 : lam[f2]);

	return p1==0.0||p2==0.0 ? 0.0 : (gamma*p1-p2)/(gamma-1.0);
}
/* compute Melbourne-Wubbena combination ---------------------------------------------
* argv   : int  fres    frequency combination
*                       (1:L12,2:L15;3:L25)
* --------------------------------------------------------------------------------- */
double Mel_Wub(const int fres,const obsd_t *rov,const obsd_t *bas,const double *lam){
	int f1,f2;
	switch (fres){
		case 1: f1=0; f2=1; break;
		case 2: f1=0; f2=2; break;
		case 0: f1=1; f2=2; break;
		default: f1=0; f2=1;
	}
	if (rov->L[f1]==0.0||rov->L[f2]==0.0||rov->P[f1]==0.0||rov->P[f2]==0.0||
		bas&&(bas->L[f1]==0.0||bas->L[f2]==0.0||bas->P[f1]==0.0||bas->P[f2]==0.0))
		return 0.0;
	double f1_f2=lam[f2]/lam[f1];
	double L1,L2,P1,P2;
	if (bas){
		L1=(rov->L[f1]-bas->L[f1])*lam[f1]; L2=(rov->L[f2]-bas->L[f2])*lam[f2];
		P1=rov->P[f1]-bas->P[f1];           P2=rov->P[f2]-bas->P[f2];
	}
	else {
		L1=rov->L[f1]*lam[f1]; L2=rov->L[f2]*lam[f2];
		P1=rov->P[f1];         P2=rov->P[f2];
	}
	return (1.0/lam[f1]-1.0/lam[f2])*
		(L1/(1.0-1.0/f1_f2)-L2/(f1_f2-1.0)-P1/(1.0+1.0/f1_f2)-P2/(f1_f2+1.0));
}
/* compute narrow-lane ambiguity -------------------------------------------------- */
double Narrow(const int fres,const obsd_t *rov,const obsd_t *bas,const double *lam) {
	int f1,f2;
	switch (fres){
	case 1: f1=0; f2=1; break;
	case 2: f1=0; f2=2; break;
	case 0: f1=1; f2=2; break;
	default: f1=0; f2=1;
	}
	if (rov->L[f1]==0.0||rov->L[f2]==0.0||rov->P[f1]==0.0||rov->P[f2]==0.0||
		bas&&(bas->L[f1]==0.0||bas->L[f2]==0.0||bas->P[f1]==0.0||bas->P[f2]==0.0))
		return 0.0;
	double f1_f2=lam[f2]/lam[f1];
	double L1,L2,P1,P2;
	if (bas){
		L1=(rov->L[f1]-bas->L[f1])*lam[f1]; L2=(rov->L[f2]-bas->L[f2])*lam[f2];
		P1=rov->P[f1]-bas->P[f1];           P2=rov->P[f2]-bas->P[f2];
	}
	else {
		L1=rov->L[f1]*lam[f1]; L2=rov->L[f2]*lam[f2];
		P1=rov->P[f1];         P2=rov->P[f2];
	}
	return (1.0/lam[f1]+1.0/lam[f2])*
		(L1/(1.0+1.0/f1_f2)+L2/(f1_f2+1.0)-P1/(1.0-1.0/f1_f2)+P2/(f1_f2-1.0));
}
/* compute ambiguity combination -------------------------------------------------- */
double amb_cmb(const double Lr,const double P){
	return Lr==0.0||P==0.0? 0.0 : Lr-P;
}
/* compute time-cross difference -------------------------------------------------- */
double single_time(const double obs_t1,const double obs_t2){
	return obs_t1==0||obs_t2==0? 0.0 : obs_t1-obs_t2;
}

/* execute command ---------------------------------------------------------------- */
int execcmd(const string StrCmd){
#ifdef WIN32
	PROCESS_INFORMATION info;
	STARTUPINFO si={ 0 };
	DWORD dwStat;
	string cmds;

	si.cb=sizeof(si);
	cmds="cmd /c "+StrCmd;
	if (!CreateProcess(NULL,(LPTSTR)cmds.c_str(),NULL,NULL,FALSE,CREATE_NO_WINDOW,NULL,
		NULL,&si,&info)) return -1;
	WaitForSingleObject(info.hProcess,INFINITE);
	if (!GetExitCodeProcess(info.hProcess,&dwStat)) dwStat=-1;
	CloseHandle(info.hProcess);
	CloseHandle(info.hThread);
	return (int)dwStat;
#else

	return system(StrCmd.c_str());
#endif
}
/* create directory --------------------------------------------------------------- */
void createdir(const string StrPath){
	string strBuff;
	size_t strPos;

	strBuff=StrPath;
	if ((strPos=strBuff.find(FILEPATHSEP))==string::npos) return;
	strBuff=strBuff.substr(0,strPos);

#ifdef WIN32
	/* convert of strBuff */
	CString CSpath=CString(strBuff.c_str());
	CreateDirectory(CSpath,NULL);
#else
	mkdir(strBuff.c_str(),0777);
#endif
}
/* uncompress file -------------------------------------------------------------------
* uncompress (uncompress/unzip/uncompact hatanaka-compression/tar) file
* args   : string  SrcFile     I   input file
*          string  UncFile  O   uncompressed file
* return : status (-1:error,0:not compressed file,1:uncompress completed)
* note   : creates uncompressed file in tempolary directory
*          gzip and crx2rnx commands have to be installed in commands path
*---------------------------------------------------------------------------------- */
int rtk_uncompress(const string SrcFile,string UncFile){
	int nStat=0;
	size_t strPos;
	string strCmd,strTmpfile,strBuff,strDir,strFname;

	strTmpfile=SrcFile;
	if ((strPos=strTmpfile.find('.'))==string::npos) return 0;

	/* uncompress by gzip */
	if (!strTmpfile.substr(strPos).compare(".z")||!strTmpfile.substr(strPos).compare(".Z")||
		!strTmpfile.substr(strPos).compare(".gz")||!strTmpfile.substr(strPos).compare(".GZ")||
		!strTmpfile.substr(strPos).compare(".zip")||!strTmpfile.substr(strPos).compare(".ZIP")) {

		UncFile=strTmpfile.substr(0,strPos);
		strCmd="gzip -f -d -c \""+strTmpfile+"\" > \""+UncFile+"\"";

		if (execcmd(strCmd)) {
			remove(UncFile.c_str());
			return -1;
		}
		strTmpfile=UncFile;
		nStat=1;
	}
	/* extract tar file */
	if ((strPos=strTmpfile.find('.'))!=string::npos&&!strTmpfile.substr(strPos).compare(".tar")){

		UncFile=strTmpfile.substr(0,strPos);
		strBuff=strTmpfile;
		strFname=strBuff;
#ifdef WIN32
		if ((strPos=strBuff.find('\\'))!=string::npos) {
			strDir=strBuff.substr(0,strPos); strFname=strBuff.substr(strPos+1);
		}
		strCmd="set PATH=%%CD%%;%%PATH%% & cd /D \""+strDir+"\" & tar -xf \""+strFname+"\"";
#else
		if ((strPos=strBuff.find('/'))!=string::npos) {
			strDir=strBuff.substr(0,strPos); strFname=strBuff.substr(strPos+1);
		}
		strCmd="tar -C \""+strDir+"\" -xf \""+strTmpfile+"\"";
#endif
		if (execcmd(strCmd)) {
			if (nStat) remove(strTmpfile.c_str());
			return -1;
		}
		if (nStat) remove(strTmpfile.c_str());
		nStat=1;
	}
	/* extract hatanaka-compressed file by cnx2rnx */
	else if ((strPos=strTmpfile.find('.'))!=string::npos&&strTmpfile.length()-strPos>3&&
		strTmpfile[strPos+3]=='d'||strTmpfile[strPos+3]=='D') {

		UncFile=strTmpfile;
		UncFile[strPos+3]=strTmpfile[strPos+3]=='D' ? 'O' : 'o';
		strCmd="crx2rnx < \""+strTmpfile+"\" > \""+UncFile+"\"";

		if (execcmd(strCmd)) {
			remove(UncFile.c_str());
			if (nStat) remove(strTmpfile.c_str());
			return -1;
		}
		if (nStat) remove(strTmpfile.c_str());
		nStat=1;
	}
	return nStat;
}